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ABSTRACT 

Variational methods (VM) sensitivity analysis employed to derive the costate (adjoint) 
equations, the transversality conditions, and the functional sensitivity derivatives. In the derivation of 
the sensitivity equations, the variational methods use the generalized calculus of variations, in which 
the variable boundary is considered as the design function. The converged solution of the state 
equations together with the converged solution of the costate equations are integrated along the 
domain boundary to uniquely determine the functional sensitivity derivatives with respect to the 
design function. 

The application of the variational methods to aerodynamic shape optimization problems is 
demonstrated for internal flow problems at supersonic Mach number range. The study shows, that 
while maintaining the accuracy of the functional sensitivity derivatives within the reasonable range 
for engineering prediction purposes, the variational methods show a substantial gain in 
computational efficiency, i.e., computer time and memory, when compared with the finite difference 
sensitivity analysis. 
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1. INTRODUCTION 


1.1 Overview of Aerodynamic Design Optimization and Sensitivity Analysis 

In the early times of flight, improvement of vehicle performance was mostly based first 
on intuition, empirically accumulated databases, and cut-and-try procedures [1,2]. Even recently, 
wind tunnel testing is being employed to perform optimization work to obtain airfoil 
performance criteria[3]. While this approach gave many valuable technical assistances, it was 
unable to furnish quick and reliable information to perform on-line design changes. 

In recent years, aerodynamic performance has been analyzed by a method of 
mathematical optimization. Eventhough there are many ways of optimization, we concenatrate 
only on the methods of optimization that require gradient information. 

1.1.1 Gradient-Based Methods (Numerical Design Optimization) 

With the advantage of modem hardware and software computer technologies, numerical 
design optimization and sensitivity analysis are currently being used to investigate the complete 
aircraft design problem using two-dimensional Navier-Stokes and three-dimensional Euler 
equations. These techniques are discussed here briefly. 

1 . 1 . 1.1 Finite Difference Sensitivity Analysis The simplest, but the most expensive, sensitivity 
analysis technique used by gradient-based optimization methods is the finite difference approach. 
This method uses the one-sided or central-difference alternative to evaluate the sensitivities of 
performance functionals, and consequently, the computational time invested would increase with 
the increment of the number of design variables. This is due to the requirement to perturb each 
design variable by an appropriate step size and then compute the flow field variable for each new 
perturbed design variable with the chosen flow solver. This approach has an additional problem 
to determine the correct step size a priori so that the correct gradient is predicted within a given 
degree of accuracy. Despite its shortcomings, Huddelston and Mastin [4], and others, have 
applied this approach in their design procedure with Euler and Navier-Stokes equations. In the 
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optimization package for general purposes optimizations, Vanderplaats [5] has also incorporated 
finite difference as an alternative to acquire the gradient information. 

1.1. 1.2 Discrete Sensitivity Analysis The other category of sensitivity analysis technique is the 
discrete analysis approach. The computation of the sensitivity equations is based on the Implicit 
Function Theorem. Due to the implicit dependence of the functional (objective function) and 
constraints on the flow field quantities, the determination of the sensitivity derivatives is related 
to obtaining the derivatives of the flow field vector with respect to the design variables. As the 
flow field equations are in most cases solved in a computational domain, the functional 
dependence of the metric terms and the coordinate points with respect to the design variables are 
also required. This approach first calls upon the multiplication and assembly of various terms to 
a very large sparse linear algebraic equations, which depend on the number of design variables, 
and then solution of these sparse system of algebraic equations for the derivatives of the solution 
vector with respect to the design variables. Despite the large computational intensity and huge 
memory requirements of this approach, the versatility to incorporate many types of constraints, 
the need to perform multidisciplinary designs of moderate geometrical complexity, and the 
flexibility to incorporate it with any existing optimization algorithm make it attractive to perform 
design and shape optimizations. 

A wealth of literature can be found for this category. Hicks [6] and Vanderplaats [7,8] 
have used the discrete approach to design airfoils in transonic flow regimes. Pittman [9] has also 
used this procedure for supersonic flow conditions. Using the small perturbation equations in two 
dimensional flows, Elbana and Carlson [10] have also employed the technique. Eleshaky [11] 
used this method for both internal and external flow problems. He also integrated the domain 
decomposition method in solving the sensitivity equations. Burgreen [12] further extended the 
methodology to the three-dimensional wing optimization and introduced an efficient way of 
parameterizing the curves and surfaces using the Bezier polynomials. With a variant of 
approximation to the fluid flow, Taylor et al. [13], Newman et al. [14], and Hou et al. [15] 
introduced an incremental iterative technique to obtain the gradient information. In doing so, 
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they have applied this new approach not only to the two-dimensional Euler and thin-layer 
Navier-Stokes turbulent equations for internal and external flows but also to the three- 
dimensional Euler equations in supersonic flow regimes. 

f .1.1.3 Variational Sensitivity Analysis The new emerging sensitivity analysis technique for 
gradient-based optimization methodology within the optimization community is the continuous 
sensitivity (variational sensitivity) analysis which fully exploits the variational methods. From 
the modified functional, this approach derives a set of partial differential equations (PDEs), i.e. 
the costate equations with their boundary conditions and the sensitivity equations. In computing 
the sensitivity derivatives with respect to the control points or design variables, this approach 
makes use of the converged solution of the state and costate equations. 

In recent years, variational sensitivity analysis has significantly contributed to the 
progress of aerodynamic design optimization. Pironneau [16] showed the usefulness of the 
variational approach in fluid mechanical problems by illustrating how to compute the minimum 
drag profile in two-dimensional viscous and laminar flows. Chen and Seinfeld [17) developed a 
methodology to compute the performance sensitivity derivatives using optimal control theory. 
Koda et al. [18] used this procedure to solve atmospheric diffusion problems. Koda [19 - 21] 
further developed this approach and outlined a numerical algorithm for the computation of 
functional derivatives. This approach is well suited to solving the optimum design problems in 
fluid mechanics. Meric [22,23] treated optimal control problems governed by parabolic and 
elliptic partial differential equations and solved them numerically using variational methods. In 
their effort to compare the gradients obtained by "implicit" and "variational" approaches, Shubin 
and Frank [24] implemented VM to optimize the shape of a nozzle of a variable cross - sectional 
area for steady one-dimensional Euler equations. Jameson [25] regarded the boundary of the 
flow domain as a control parameter and then designed airfoils using the potential as well as the 
two- and three-dimensional compressible inviscid flows. Cabuk and Modi [26] implemented a 
perturbation method to compute the optimum profile of a diffuser for a maximum static pressure 
in a two-dimensional steady viscous incompressible flow. Ta’asan et al. [27] have successfully 
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implemented variational methods and optimized an airfoil in the potential flow field. Quite 
recently, Ibrahim and Baysal [28] demonstrated the versatility of the variational methods to solve 
aerodynamical design problems for internal flows in different Mach number regimes including 
shock flows. Following the same approach as Jamson [25], Reuter and Jameson [29] optimized 
airfoils in potential flows. Iollo and Salas [30] used variational methods to solve two- 
dimensional internal flow optimization problem with embedded shock to match a pressure 
distribution. In this class of optimization, the functional sensitivity derivatives are directly 
coupled to the solution of a set of linear partial differential equations, i.e., the costate equations 
and their boundary or transversality conditions that result from the variation of the augmented 
Lagrangian function. The success of any optimization by this approach is, therefore, destined to a 
stable and converged solution of the costate equations. 

1.2 Generality of the Variational Approach 

First, since the costate equations are once and for all derived from the continuous PDE of 
the state equation, any robust solution method can be adopted to furnish the converged solution 
so that the costate equations can be solved until convergence is attained. This means that one 
does not necessarily have to solve the original state equation from which the costate equations 
are derived. Secondly, any other convenient discretization methods different from the type of 
discretization one uses for the state equations can be selected for the costate equations. The 
requirement that the costate equations be discretized exactly the same way as the state equations 
is shown not to be necessary, at least for quasi one-dimensional Euler equations [31]. Thirdly, 
any time integration method different from the time integration method used for the state 
equation can be selected to advance the costate equations to steady state. The fourth point to 
mention is the design variables. In the approach proposed, note that the shape of the domain is 
considered as the design parameter, and its contribution to the functional sensitivity derivatives is 
directly incorporated as shown in Ref. 32. 
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2. AERODYNAMIC DESIGN OPTIMIZATION AND SENSITIVITY ANALYSIS 

2.1 Constrained Optimization Methodology 

A constrained optimization method in general encompasses three elements of 
optimization, i.e., design variables, constraints, and objective function. These are discussed here 
briefly. 

2.1.1 Design Variables in Variational Sense 

In most aerodynamic optimization problems, the design variables are generally of a 
geometric nature, such as the coefficients of some geometric functions, surface grid points [11], 
aerofunctions [33], or polynomial functions such as Bezier-Bemstein functions [12,34] and 
spline functions [35]. 

Variational methods treat the boundary of the domain in a continuous fashion, and 
therefore, the boundary is considered as part of the solution to the design problem. With the 
assumption that the domain Q is sufficiently regular, the location of points on the boundary X v 
can be considered as a continuous design variable. Mathematically, the coordinates of the 
varying boundary in the continuous sense can be expressed as 


X r =/(*„) (>) 

where X d are the design variables. In aerodynamic optimization problems, the vector of design 
variables is provided for very limited and simplified geometries, for instance, four digit NACA 
airfoils and some nozzles. However, for general-purpose geometries, these control points must 
be determined through iterative methods from certain functional relationships such as the Bezier- 
Bemstein polynomials [12]. Because these polynomial functions are known to generate smooth 
curves and surfaces for a minimal number of control points, the function / which describes the 

curve for the two-dimensional problem, is given by [34] 
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f(u) = lX dJ B hn (u) 

i-O 

where 


for u G:[0,l] (2) 


B i n (u) = C(n,i)u‘(l - t7)"" 


(3) 




n\ 

/!(«-/)! 


(4) 


In Eqs. (2) - (4), B in (u) are the blending functions, which are key to the behavior of the curve, 
C(n,i) are the binomial coefficients, u is the normalized arc length and n is the order of the 
Bezier-Bemsten polynomials. In this study, with the use of Eqs. (1) - (4), the location of the 
control points can be considred as the design variables. 

2.1.2 Constraints 

Constraints are the integral parts of the optimization procedure that influence the final 
outcome of the functional. They can be geometrical, flow-type, equality or inequality constraints, 
or a combination of all or some that depends on the particular optimization problem one wants to 
address. 

In the variational formulation of design optimization problems, the flow-type constraints 
are expressed in the integral forms. The geometrical and side constraints, on the other hand, can 
be formulated either in the integral or discrete forms. For the general variational approach, 
generic flow-type constraints are expressed as 

Gj(Q,n) = [gj(Q,n)dT <; 0 fory - 1, 2, ..., nconf (5) 

r 


where T , is the deformed boundary and nconf is the number of generic fluid-type constraints. 
The generic geometric-type and the side constraints can also be given as 
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Gj(X d ) s 0 for j = nconf+1, nconf+2, ..., neon (6) 

and 

Xj™ < Xjp < ^7" for i = 1 , 2, ndv (7) 

where neon is the total number of constraints, and ndv is the number of design variables, 
respectively. 

2.1.3 Objective Functional 

In the variational methods (VM), the objective functional is defined in the form of a 
definite integral involving an unknown state function Q, which can be dependent on some 
normal vectors n and other problem parameters. Then, the objective functional is extremized at 
the converged state solution over the curve of the surface described by the vector of design 
variables. Mathematically, a generic functional on the boundary J f , is defined as 

J T (Q,n)-[D(Q,n)dT (8) 

r 

where D, for the two-dimensional problem, is the objective function specified on the curve or 
boundary. The selection of the objective function is mostly dictated by the flow physics. 

2.2 Variational Formulation of Aerodynamic Optimization Problem 

When constraints are involved in the optimization problem, the partial derivatives of the 
functional and the constraints cannot be zero at the same time since they are functionally related 
to each other through the optimality criteria [36, 37], One common practice is to cast the 
constrained optimization to unconstrained optimization through the introduction of the weighting 
functions or Lagrange multipliers X(X). The other is to sequentially solve a linear or quadratic 
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programming problem, which is an approximation of the original constained minimization 
problem. In the later approach, one needs to derive the sensitivities of the performance 
functional. 

To start the derivation, the steady state solution of the two-dimensional Euler equations, 
i.e., the residual R(Q), is written as 


R(Q) = E X+ F y = 0 


( 9 ) 


and the generic boundary conditions are expressed as 

H(Q,n) = 0 (10) 

Without changing its value, the objective functional J f can now be modified as 

J fa (Q,n,X,fi) = J f +[X T (X)R(Q)dQ + ^ T H{Q,n)dT (11) 

q r 

where T and Q are the deformed boundary and domain, and A and j2 are vectors to be 
determined. 

2.2.1 Standard Formulation of an Aerodynamic Optimization Problem 

A mathematical formulation of the constrained optimization problem can be expressed as 

min |j r } (12) 

subject to 

Gj(Q,n) = f gj (Q,n)dr ^ 0 


for j =1,2, ..., nconf 


(13) 
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Gj{X d )<. 0 j = nconf+1, nconf+2, neon (14) 

and 

X‘°J er <; X ID z X%’ er i =1,2, NDV (15) 

where the flow field variables Q are the solution to the state equations, R(Q) ■ 

2 . 2.2 Derivation of Functional Sensitivity Equations 

As the detailed derivation for functional sensitivity is showen in Ref. 32, we give here 
only the high lights of the procedure. Hence, the variation of the fluxes (to the first order) can be 
written as 

6E = A6Q and 6F = B6Q (16) 

where A , and B are the Jacobian matrices in the x and y directions, respectively. Then the 
fluxes on the deformed space due to the variation of the boundary can be approximated as 

E(S)-E(Q) + SE(Q) and F(§) = F(0 + 6F(0 (17) 

By application of the principles of calculus of variations, the variation of the modified functional 
can be approximated by [36] 

A J T „ - ),?(£))- (18) 

where X and X are position vectors of the deformed and undeformed coordinate systems, 
respectively. Then, the Taylor expansion of the integrand of Eq. (18) is computed (the linear part 
relative to e ) as [36-38] 
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5J Va = jjl) + D 0 8Q + D n 6n + DKdn^dT + f(fi r + <5/2 r ) j// + H q 8Q + H n 8n + HKdn\dT 

+ j|(A r + 6xf[E x (Q) + (8E(Q)l + ^(0 + (<5F(0),]}|J S |</Q 
- jA r [l,(0 + F y (Q)]dQ-^D + /Ttfjr (19) 


where £ is a small parameter, J s is the space transformations matrix that is given as 
/ + V. [32] and the quantity k in Eq. (19) is the curvature and can be calculated as [37] 


K = -V°n (20) 

where ° denotess a dot or inner product and n is the unit normal which can be computed from 
the grid generating routine or from the analytical derivatives of the Bezier-Bemstein polynomials 
as [34] 




3 ' 

du du~ 

du 



where ® is a vector multiplication sign and f(u ) is defined in Eq. (2). In Eq. (19), 5n is defined 
as 

6n = 6Xon ( 22 ) 

Now, by taking only the linear terms of Eq. (19), one obtains 

8J ra = jjz) 0 <5(? + D n 8n + Dk5h + /2 7 |// ( . ; b(9 + H n 6n + //x<5«jj<7T 

+ j{l r [(<5£(0), +(.SFm y }d& + /{(^n^<2) + F^da+fS^HdT 
+ JA T \I X + F.JSnrfT (23) 
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Using Eq.(16) and performing integration by parts, the second term in Eq. (23) is expressed as 


j{i r [(«s(0), +(«F(fi)),}<M=/{-&4 -i r y B\5Qda 


+ 


l{^[ An - 


+ Bn v ]dQ)dT 


(24) 


Substitution of Eq. (24) into Eq. (23) gives 


&J, 


r a 


-•f{ 


D 0 dQ + D n 5n + Dk6h + fl 


^ t [h q 6Q 


+ H„dn + HkSh 


■dr 


+ 


j{a T [An x + Bn y ]dQ \dT + j{-A^ - A^}<5g</Q 
r ^ J q 

jf|(6A) r [^(0 + F y {Q)^d£l+ p^Hdr + JfA r [E x + F^ndT (25) 


Note that for the arbitrary variations of <5A and 6/x and with Eqs. (9) and (10), the last two terms 
in Eq. (25) are identically zero. Then Eq. (25) reduces to 

6J ro = j| D q 6Q + D„6n + Dk5h + fl T \ji Q d§ + H n 5n + 77x<5n|tfT 

+ ffr[An x + Bn y ]dQ]dT+ ^\-X t x A - k T y B^dt)dQ+ fr T [l x + F^ndT 
r ^ J q 1 r 

(26) 


In Eq. (26), the vectors A , and jx can now be determined to eliminate the terms 
associated with 8Q . Consequently, the costate (adjoint) equations are given as 

-A r X - B r k = 0 in Q (27) 

y 


Upon the combination of Eqs. (26) and (27), the variation of the functional becomes 
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6J ra - jj. D q 5Q + D n 8n + Dk6h + fi T ^H Q 8Q + H„6n + Hicdn 

+ $\*. T [An x + Bn y ]d§\dT + J*A 7 '[X + F^ndT 
xi'- * r 


■dr 


(28) 


With Eq. (A.l 1) in Appendix A of Ref. 32, we now express (5 <2 in terms of 6Q to get 


5Q=8Q-^dX 

dX 


(29) 


For the sake of computational simplification, the variation of dX on the boundary is limited only 
to the y component in this study, i.e.. 


8X = [0,dy]' 


(30) 


and Eq. (29) is simplified as 


SQ.SQ-^-dy 

dy 


(31) 


Also an approximatation of Eq. (22) and use of Eq. (30), Eq. (22) can be written as 


6n = 6X o n 


= [0,<5j>f °[n x ,« y ] 


= n y °dy 


(32) 
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By use of Eqs. (29) - (32), Eq. (28) is now given as 

SJ ra = J \d q 5Q + n y D„8y + n y DKdy - n y D Q Q y 6y +}</T 
r 

+ J ^j2 T \jI Q dQ + n y H„dy + n y Hicdy - n y H Q Q y 8y^dT 
r 1 

.+ /U r [In, + Bn y \dQ - n y Q y 6y)\dr + + F^SydT (33) 

r 


For arbitrary 5Q and the variation of y on the boundary T, Eq. (33) gives the boundary 
conditions for the costate equations and the sensitivity equations, respectively, as 


r— V r 


An r + 


B /Tyj X + Dq + 


H T Q ji\ - 0 


on r 


(34) 


and 


6J ra = j| D„ + Dk- D Q Q y + fl r ^H n + Hk - H Q Q^n y dydT 
r v 

- (\i T Un x + S*,]e, \n y 6ydT + f A r [E, + F y }, y dydr (35) 

r 1 J r 

The unique determination of Eq. (35), therefore, demands the unique and converged solutions of 
Eqs.(9) and (10), (27), and (34). 

2.2.3 Derivation of Constraint Sensitivity Equations 

With the constraints defined in Eq. (13), the residual, and the boundary conditions, Eqs. 
(9) and (10), one can formulate the modified constraints as 

Oja**,,*,) - fg,(a»yr + U/[E X + + lH/H(Q,n)dT 

r Q r 


for j =1,2, ..., nconf 


(36) 
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By following the same procedure as was done for the objective functional, the costate equations, 
boundary conditions, and the constraint derivative coefficients, respectively, can be expressed as 


| -A T X jx - = 0 for j = 1, 2, ..., nconf in Q 

| \An x + Bn y Jkj + g] Q + Hgiij J = 0 for j = 1, 2, ..., nconf on T 

6G Jr a =/{^„ + Sj k - SjoQy + Vj T [H n +Hk- H Q Q^n y SydT 
r 

- r{X/[I«, + S„,]e r Lsydr + fxp, + F y \ r 6ydr 

for j =1,2, ..., nconf 


(37) 

(38) 


(39) 


As can be discerned from Eq. (39), the computation of the constraint sensitivity equations 
requires the solution of a new set of costate equations and boundary conditions as many times as 
the number of constraints. 


2.3 Numerical Optimization 

Two steps are essentially followed in this approach. The first step is to determine the 
search direction, S, and the second is to compute the magnitude of the step size a. These two 
quantities can be computed as proposed in Refs. 39 and 40 A typical computation of the feasible 
direction starts at the boundary of the feasible domain, and its magnitude and directions are kept 
constant as long as the search direction keeps the design variables in the feasible domain while 
improving the performance index. Otherwise, a new search direction and step size are 
recomputed with the new gradient information and this process continues until the optimality is 
met. Mathematically, the feasible direction can be formulated as 

S T oVg.<0, 


(40) 



15 


where i is part of the active constraints and the usable direction at a point is given by 
S T oVJ r zO (41) 

The change in design must be sought along the combination of the useable and feasible 
directions so that the functional or the performance index is reduced as much as possible, and the 
design is kept away from the constrained boundary as much as possible. By the use of Eqs. (40) 
and (41) in the method of feasible directions, the new design variables are updated as 

X ^ + 1 = X” D + aS (42) 

where n is the iteration number. The values of the design variables are continuously altered until 
the criteria for the optimal solution of the performance index are satisfied. 

3. COSTATE EQUATIONS AND SOLUTION METHODS 

3.1 Introduction to the Numerical Integration of Costate Equations 

The coefficients of the costate equations are constant matrices whose components are 
derived from the converged solution of the state equations. They are globally constant in time 
and locally constant in space. But the interpretation of constant matrices must be understood in a 
sense that, during the time integration of the costate equations, only the costate variables evolve 
in space and time to convergence. The costate equations are identical to the Euler equations in 
form, but mathematically, they are different in the sense that they do not meet the homogeneity 
requirement to put them in a conservative form like the Euler equations. From the numerical 
view point, one can adopt any solution algorithm, which is used for the Euler equations, to the 
costate equations. This can be explained by the fact that the fluxes on the cell faces or at grid 
points can be artificially constructed by approximating the solution vector of the costate 
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equations either on the cell face from the right and left sides of the cell centers or at the grid 
points in exactly the same way one does for the state fluxes and solution vector. 

The costate equations, like the state equations, are solved by use of the time dependent 
techniques. The Eqs. (27) and (37) are, therefore, modified to include the unsteady term with the 
proper signs so that this time dependent technique is fully exploited. Thus, for instance, Eq. (27) 
in the generalized coordinates system is expressed as 

±A T - A r X s - B r A„ = 0 (43) 


The proper sign selection of the time term is dependent on the complementary property of the 
well-posed boundary conditions of the state and costate equations. For Eq. (43) to be well-posed, 
the positive sign of the time term is selected, and Eq. (43) becomes 

A t -i r A r 5 r A,=0 (44) 


3.2 Boundary (Transversality) Conditions 

The objective functional boundary conditions, i.e. Eq. (34), in their general forms are 
again for the sake of convenience presented here as 


| A n x + Bn v J A + Dq + i/Jjuj = 0 


on T 


(45) 


The objective functional and the no-mass penetration conditions are defined only on the solid 
boundary, and hence their derivative contributions in Eq. (45) are identically zero. Therefore, the 
boundary conditions for the inlet, exit, and center-line reduce to 


<T T F 1 

V A inlet ’ center ’ exit/ 


(46) 


+ 5« v ] r Aj = 0 


on 
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For the supersonic flow, the inlet condition is known, and hence the variation of the 
vector of the flow field is identically zero. Therefore, with Eq. (46), the values of the costate 
variables at the inlet boundary can be approximated from the internal stencils. Because the vector 
of the flow field is computed from the internal grids at the exit plane, Eq. (46) gives four linear 
independent equations for the costate variables, which result in all the costate varibles to be 
identically zero. On the centerline, the normal velocity is known to be zero, and one of the 
costate variables, for instance A 3 , is assigned a value, and the remaining flow field quantities are 

to be detemined from the resulting 3x3 system of equations as given in Eq. (45) 

One way of treating the boundary conditions, i.e. Eq. (45), is to use Eq. (10) and to find a 
relationship between the conservative field variables Q by taking the variation of Eq. (10). This 
procedure eliminates the constant Lagrange multipliers u and modifies the functional sensitivity 
derivatives, Eqs. (35) and (39) by a term resulting from the variation of the normal vector n at 
the solid boundary. 

On the solid boundary, on the other hand, the costate variables are determined by use of 
the complete form of the compatibility relationships and the sign of the eigenvalues of the costate 
Jacobian matrices. Once the values of the costate variables on the solid boundary are computed, 
the constant Lagrange multipliers j2 of the no-mass penetration condition can be calculated by 
solving the complete set of the boundary condition. The results presented in this study are 
obtained by solution of the complete boundary conditions as given in Eqs. (34) and (38). 

3.3 Linearization of Costate (Adjoint) Equations 

By the same linearization procedure we used for the state equations, Eq. (44) can be 


approximated as follows: 


k - Az ~{ 2 ’k + &k} ' {^k + 

(47) 

k - Ar^ |><5 s + B' A,]X} - {A T X t + B r A,}” 

(48) 
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A - At— jj 
T <?A U 


A t 6 § + B t 6 v 


B£-( 


A ™ — A — 

^ r A ? +5 r A^ 


(49) 


By approximation of the time and space terms, Eq. (49) becomes 

AA" - At [2 t 6^ + B r S^ AA" = At^^ + B r 6^ a}" (50) 

3.4 Time-Integration Method 

In this study we have used the implicit, i.e., the ADI method to drive costate equations to 
steady state. For the implicit method, the ADI factorization of Eq. (50) is used to split it into the 
§ and rj sweeps. Let us define the right side R x of Eq. (50) as 


R, 


- Ar{[(i-) r a E * +(>) r <5 5 - +(B-) T s; 


(51) 


where R x is the residual for the costate equations. Also, Eq. (50) can be put in its split form of 
Jacobians and fluxes as 


7- Ar[(i + ) 7 'd- +{A~) T dl +(B + ) t 6~ + (7T) 7 ' 




AA" =7?, 


(52) 


Then the £ and t) sweeps of Eq. (52) are given as 


{/-Ar[(i‘) r d-+U-) r «;]}AF -R x 

and 

{ 7 - Ar[(5 + ) r <5; + ( 5 -) 7 ' 5 ;]}aA" = AA"* 


(53) 


(54) 
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4. RESULTS AND DISCUSSION ON DESIGN OPTIMIZATION OF INTERNAL 
FLOWS USING TWO-DIMENSIONAL EULER EQUATIONS 

The main thrust of this section is to briefly discuss the the numerical results of the 
variational sensitivity analysis that are obtained by the use of two-dimensional Euler flow 
equations. Additionally, the efficiency and accuracy of the variational sensitivity in comparison 
to the finite difference are analyzed. 

4.1 Two-Dimensional Nozzle Optimization Problem Formulation 

At least a couple of reasons can be given for choosing the two-dimensional nozzle 
geometry in order to demonstrate our point of optimization methodology. The first one is that 
one can easily obtain various types of nozzle geometries by simply using already known 
analytical expressions for different flow conditions. The second important reason is also the need 
to develop a scramjet nozzle afterbody for the High-Speed Civil Transport. The third one is the 
need to develop efficient wind tunnels with optimal shapes for various experimental wind tunnel 
applications. The optimization problem demonstrated here seeks the optimal shape for the 
maximum thrust in conjunction with the nonreverse flow condition at the exit. Hence, the 
example problem is formulated as the maximization of the functional defined by 


J f =[PdT (55) 

r 

with the constraint that the static pressure P at the exit assumes a certain percentage of the 
ambient pressure p x for maximum expansion at that exit lip of the solid boundary. Therefore, 

the constraint is mathematically posed as 



k/T<;0 


(56) 
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4.2 Two-Dimensional Nozzle Flow 

The initial geometry for this internal flow configuration is given in Fig. 1. It is a 
supersonic nozzle where only half of the physical domain is considered with 137 x 69 grid 
points. It is a convex type of geometry with the smallest area at the inlet and a diverging 
afterbody for supersonic expansion. The only aerodynamic inequality constraint considered is the 
criteria on the static pressure at the exit lip of the nozzle to reach a certain percentage of the 
ambient pressure as a necessary condition to avoid any reverse flow from underexpansion as the 
shape evolves during the optimization cycle. 

To assess the variational methods for sensitivity analysis, computational efficiency and 
accuracy calculated by variational methods and finite difference are compared. One of the 
obvious limitations with the finite difference is the uncertainty to a priori determine the step size 
that will give reliable sensitivity derivatives. The magnitude of the stepsize is dependent on how 
accurate one needs the derivatives to be. If, for instance, one only needs a 10% deviation from 
the assumed exact derivative, then the step size must be under a 10% range of the derivative. In 
our case of computing the sensitivity derivatives using the finite difference, we have assigned the 
step size to be 0.0001. 

The x component of the design variables (Bezier control points) are a priori computed 
as being spatially invariable, and the variation of the design variables is allowed only in the y 
direction. This apparent limitation of the design variables must not be a hindrance since addition 
or deletion of any desired design variables in the design domain will produce the same result. To 
verify this claim, two sets of design variables, in addition to the assumed optimal number of 
design variables (in this case the optimum is eight design variables), were investigated. The first 
set was performed by increasing the number of design variables by four and the second one by 
decreasing it by four from the optimal number of design variables. Here, the optimal number 
defined as that number of design variables which reproduces the closest shape to that of the 
initial geometry. 
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As presented in Table 1, the CPU time and memory requirements of complete cycle of 
optimization for the two additional sets of design variables are almost identical for the two- 
dimensional optimization case. Therefore, the eight design variables are considered as the 
optimal number of design variables which produced the desired computational efficiency for our 
test case. On the other hand, this slight memory increase as the number of design variables 
increases could be a warning to the eventual computational memory increase as the 
dimensionality, number of constraints, and design variables increase. The second aspect of the 
role played by the number of design variables may be the influence on the optimal shape (Fig. 2). 
All three categories of the design variables produced slightly different optimum shapes from 
each other. Comparing all three shapes (Fig. 2), the shape produced by twelve design variables 
appears to follow the shape produced by the four design variables in the compression area 
(upstream) and the shape of the eight design variables in the expansion area (downstream). The 
shape generated by the optimal number of eight design variables shows a slight change of shape 
upstream, from approximately x=0.1 to x=0.375, and downstream, from approximately x = 0.7 to 
x = 1 .0 as compared with the shapes generated by the other sets of design variables. The shape 
change in the compression area seems to be more desirable beause it produces high-pressure 
ratios and thereby gives more thrust as one integrates the change of pressure along the changing 
nozzle shape. The shape change in the expansion region, on the other hand, reduces the ratio of 
the static pressure to the ambient pressure, which results in less thrust augmentation. This 
physical phenomena is further reflected in Fig. 3 where the optimal thrust of the eighth design 
variable shape is higher than the other two design variable shapes. 

From the parametric studies (four, eight, and twelve design variables) conducted, one 
may conclude that the eight design variables are the optimal number of design variables to 
sufficiently represent the nozzle shape and at the same time to give a better thrust and 
computational efficiency. 

The evolution of the design variables for the variational methods and the finite difference 
approach are given in Figs. 4 and 5, respectively. Except for the second and the seventh design 
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variables, the general trend of the evolution of the design variables in both approaches is similar. 
In the variational case, the second design variable approaches the first design variable and the 
seventh one tends to come close to the eighth design variable. In the finite difference case, 
however, the second and the seventh design variables tend to pull away from the first and eighth 
design variables, respectively. As shown in Fig. 6, due to the movement of the second and the 
seventh design variables in the opposite direction, the optimal shapes of the variational methods 
and finite difference are slightly different. As explained in the parametric studies, the decrease of 
the optimal (as compared with the initial) shape or optimal design variables in the compression 
region is much more advantageous to the decrease of the optimal shape or optimal design 
variables in the expansion region for the supersonic flow regime. This is due to the effect that the 
decrease of the shape in the upstream results in the substantial gain of high pressure ratio 
(compare Figs. 7 and 8) which favors the augmentation of more thrust (Fig. 9) in the design 
process. Figure 9 also clearly indicates that the pressure distribution in the expansion region in 
general and at the lip of the nozzle in particular is within the constraint specification as imposed 
in the aerodynamical constraint given by Eq. (56). 

As given in Table 2, the accuracy of the variational methods is verified by comparing the 
variational functional sensitivity derivatives to the functional derivatives of finite difference. If 
one takes into consideration that the sensitivity coefficients of the finite difference are dependent 
on the step sizes, then the gradient values obtained by the variational methods are well within the 
engineering prediction range, except for the second and the seventh sensitivity coefficients. The 
discrepancy of those two sensitivity values may be associated with the difficulty to properly 
implement the boundary conditions of the adjoint equations. Despite the differences on these two 
sensitivity derivatives which correspond to the second and seventh design variables, the optimal 
shape and thrust of the variational methods are comparable with those of the finite difference as 
presented in Table 3 and Fig. 6. It is known that the finite difference uses function evaluations to 
compute the gradient information while the variational methods solve another set of partial 
differential equations and sensitivity derivative equations. Due to this, there is a memory 
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increment of approximately 1.3 mega words as shown in Table 4. This slight increment in 
memory is negligible as compared with the other gradient-based sensitivity analysis approaches, 
such as the discrete sensitivity analysis which requires higher memory allocation for the given 
optimization problem. 


5. CONCLUSIONS 

A two-dimensional nozzle optimization problem was considered, and the application of 
variational methods to compute the optimal shape for the maximum thrust is presented. During 
the design process, the supersonic nozzle remained supersonic while improving the performance 
index or thrust (Table 2). Also, while the VM's computational accuracy (Table 3) is comparable 
with the finite difference, its computational efficiency and memory savings (Table 4) are found 
to be substantial. As memory and computational efficiency are the bottle-necks for large two- 
dimensional and three-dimensional problem in general, variational methods are one of the most 
viable candidates in solving design optimization problems. 
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Table 1. CPU Time and Memory for Four, Eight, and Twelve Design 
Variables With Variational Methods 


Design variables 

CPU time (sec) 

Memory (MGW) 

4 

868.0463 

5.249459 

8 

864.2226 

5.249939 

12 n 

866.2128 

5.250579 


Table 2. Sensitivity Derivatives by Variational Methods and Finite Difference 



Variational methods 

Finite 

difference 

Deviation (%) 

1 

9.1483E-2 

9.444 IE-2 

3.1 

2 

7.9228E-2 

1.1062E-2 

86.0 

3 

-6.6563E-2 

-4.7906E-2 

28.1 1 

4 

-5.549 IE-2 

-5.9409E-2 

6.6 

5 

-4.642 IE-2 

-5.3278E-2 

12.9 

6 

-3.8979E-2 

-4.6287E-2 

15.8 

7 

-3.2186E-2 

-6.9988E-2 

53.9 
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Table 3. Initial and Optimal Values of Functional and Constraint 
for Variational Methods and Finite Difference 



Variational methods 

Finite 

difference 

Initial 

Functional 

0.045481 

0.045481 

Constraint 

-2.10787 

-2.10787 

Optimum 

Functional 

0.049958 

0.049885 

Constraint 

-0.5858 

-0.5668 


Table 4. Efficiency Comparison Between Variational Methods and Finite Difference 



Variational 

methods 

Finite 

difference 

CPU 

time (sec) 

Complete 

optimization 

865.098 

4356.33 

Single 

analysis 

Euler equations 


128.23 

Co-state 

equations 

58.59 


Memory 

(MGW) 

Complete 

optimization 

5.25 (with 
sensitivity 
eqs.) 

3.98 (no 
sensitivity 
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